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ABSTRACT 


The theoretical development and application of 
quasilinearization techniques to problems of system 
identification are presented. A one-degree-of- 
freedom system is used for the examples for the 
identification problems. Several numerical experi- 
ments are presented that relate the effects of noisy 
input data to the accuracy obtained in computing 
unknown parameters of the governing differential 
equations. Also presented are the application of 
quasilinearization to the identification problems of 
automobile coasting and parachute dynamics. Con- 
clusions from these experiments are included. 


APPLICATIONS OF QUASILINEARIZATION THEORY 

TO SYSTEM IDENTIFICATION 

By George A. Zupp, Jr. , and S. Bart Childs* 
Manned Spacecraft Center 


SUMMARY 


The theoretical development and application of several quasilinearization tech- 
niques to system- identification problems are described. Presented in the theoretical 
development of the quasilinearization techniques are (1) the definition of system iden- 
tification, (2) the development of the Newton- Raphson method for solving nonlinear si- 
multaneous algebraic equations, and (3) the extension of the Newton- Raphson method by 
Kantorovich to the solution of nonlinear or linear differential equations subjected to 
multipoint boundary conditions. 

System- identification problems are solved for four one-degree-of-freedom sys- 
tems, a linear and a nonlinear oscillator, a free-falling parachute, and a coasting auto- 
mobile. Several numerical experiments are presented to relate the effects of noisy data 
to the accuracy with which the coefficients and initial conditions of the governing differ- 
ential equations can be computed. Experimental data were available for the free-falling 
parachute and the coasting automobile. The aerodynamic drag-area parameter was de- 
termined in both instances by quasilinearization. Numerical examples are presented 
to relate the accuracy of the predicted aerodynamic drag-area parameters to inaccura- 
cies in the experimental data for the free-falling parachute and coasting automobile. 


INTRODUCTION 


The analysis of a physical phenomenon usually starts with the postulation of an 
equivalent mathematical model, which is usually in the form of a differential equation 
(or equations). In differential- equation form, the coefficients and initial conditions 
usually have physical significance, and a knowledge of their values can lead to a better 
understanding of the physical phenomenon. Determination of the initial conditions and 
coefficients (unknown parameters) of the governing differential equations by "closed 
form" analytical methods is not always possible. In certain instances, measured ob- 
servation of the phenomenon can be used in determining the unknown parameters. De- 
termination of the unknown parameters of governing differential equations by using 
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measured observations of the phenomenon is called system identification. This report 
presents the development and application of several quasilinearization techniques to 
the solution of system-identification problems. Consideration is given to the effect of 
"random" noise on the identification of nonlinear systems. Additive noise or bias er- 
rors in the data are not considered. 

Standard system-identification techniques for dynamics problems commonly rely 
on sophisticated instrumentation to measure the state variables necessary to identify 
the system. The accuracy to which a state variable can be measured is dependent on 
the measuring technique; however, the higher the order of the state variable being 
measured, the noisier the measurements usually are. Any estimates of the unknown 
parameters, based on experimental data of the system- identification problem, reflect 
these experimental errors. Therefore, experimental data of highest accuracy should 
be used in the prediction of the unknown parameters of the system- identification prob- 
lem. 


Two quasilinearization techniques can be used to calculate the unknown param- 
eters of the system-identification problems that contain experimental data. These 
quasilinearization techniques essentially determine the coefficients and initial condi- 
tions of the governing differential equations that "best" satisfy the experimental data. 
These techniques are (1) the Newton-Raphson method for solving sets of nonlinear al- 
gebraic equations, which are the solutions to the governing differential equations, and 
(2) the Newton-Raphson-Kantorovich method for solving system-identification problems 
for which closed-form solutions of the governing differential equations are not possible. 
Ideally, either technique reduces the instrumentation requirements to selection of the 
state variable that can be measured most accurately. 

The application of the quasilinearization techniques is illustrated by identifying 
the unknown parameters of several one-degree-of -freedom systems. The examples 
presented apply to the dynamics of a linear and a nonlinear oscillator, a free-falling 
parachute, and a coasting automobile. Several numerical experiments are presented to 
relate the effects of noisy boundary conditions to the accuracy with which the coeffi- 
cients and the initial conditions of the governing differential equations can be computed. 


SYMBOLS 


A, B, C, D, a differential- equation coefficients 

C coefficient matrix 


r S 
D 


aerodynamic drag-area parameter 


d inversion vector of the coefficient matrix and a 

T 

E = C C 

— T— 
e = C d 


2 


g 


J 

L 

M 

n 


2 

acceleration of gravity, ft/sec 

Jacobian matrix 

number of boundary conditions 

vehicle mass, slugs 

index 


P 


R(x) 

R (x) 
m ' 


S 

t 

x 

x 0 , Xq, y q; y q) u, ^ , x 


number of dependent variables in Z 
residue-equation vector 

residue equation of the form R m (x) = ^ m ^ x j’ x 2 ’ 

individual boundary conditions 
time, sec 

displacement (height), ft 
state variables 



= 0 


y 

z 

z 

a 

P 


dependent condition variable 

vector of dependent variables Z 
state-variable component 

parameter specified to ensure satisfaction of the boundary con- 
ditions 

rolling friction coefficient 

3 

air density, slugs /ft 


Subscripts: 

0 initial 

i, k, m, n indices 
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particular solution 
homogeneous solution, i = 1, 2, 3 
transpose 


Superscripts: 

( 0 ) 

(i) 

T 


first derivative with respect to time 
second derivative with respect to time 


Operators: 

n 

(••) 


SYSTEM-IDENTIFICATION TECHNIQUES 


System identification is the process of fitting the governing differential equations 
to a given set of boundary conditions. The process involves the determination of the 
coefficients and initial conditions that satisfy the given set of boundary conditions. The 
coefficients and initial conditions of the differential equations are referred to as the un- 
known parameters of the mathematical model. 

To illustrate the concept of system identification, the example is considered in 
which the governing differential equations are 



(1) 


and 



( 2 ) 


where a is an unknown parameter. The solution of equation (1) can be determined by 
straightforward integration. If the range of interest is from zero to t, then the solu- 
tion to the governing differential equation can be written as 


y 



( 3 ) 
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where is the initial condition. If the desired solution must satisfy specified bound- 
ary conditions such as y = y^ at t = t^ and y = t = t g , then 


a = 




and 


y 0 = y 2 ex P 





(4) 


(5) 


Thus, system identification is straightforward for a system that can be described by a 
differential equation as simple as equation (1). Although equation (1) is linear, the 
identification problem is nonlinear. 

In a more complex differential equation such as 


^ +X 2 y = 0 (6) 

dr 


the apparent solution 


y = A sin Xt + B cos Xt 


(7) 


where A and B are arbitrary parameters, is again straightforward. If at t = 0, 
y = y Q and y = y Q , then 


A 



(8) 


and 


B = 


y 0 


(9) 
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Equation (7) can now be written as 


y = y n cos Xt H — r— sin At 


( 10 ) 


To determine y^, yg, and A, the following three independent boundary conditions for 
y(t) are necessary: 


l. y = y 1 

at 

t = 

2. y = y 2 

at 

t = 

S. y = y 3 

at 

t = 


1 

2 

3 


Use of the data points tj, t^, and t g and equation (10) gives 


y o 

y 1 = y 0 cos At x + — 

„ y o 

y 2 = y 0 cos At 2 + ~y 

„ y o 

y 3 = y 0 cos xt 3 + ir 


sin Atj 
sin At 2 
sin Atg 




> 


J 


(ID 


which are independent equations if the spacing of the data points is properly chosen. 
The presence of the sine and cosine functions indicates that this set of equations is 
nonlinear. To determine y Q , y Q , and A, an iterative technique (such as the Newton- 

Raphson method discussed in another section of this report) for solving nonlinear equa- 
tions is required. 

The governing differential equations are always written as a set of first-order 
equations, and unknown constant parameters are incorporated by adding a null equation. 
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Thus, equation (6) is written as 


dy _ z 

dt ~ z 


dz .. 
dF = -* y 


M =0 

dt 


A 


> 


y 


( 12 ) 


where 


£ 



(13) 


From the discussion of equations (1) and (6), it is obvious that the problem of 
system identification is difficult when attempted analytically. With the aid of the digital 
computer, numerical techniques are available for the solution of complicated system- 
identification problems. A discussion of quasilinearization techniques and their appli- 
cation to system identification are presented in the following sections of this report. 


QUASILINEARIZATION TECHNIQUES 


The concept of quasilinearization is presented in the theory of dynamic program- 
ing (ref. 1). Quasilinearization techniques provide a systematic, iterative approach 
for solving linear and nonlinear algebraic and differential equations. Taylor and Ilift 
(ref. 2) have successfully applied a Newton-Raphson method to the solution of the 
system- identification problem of determining aerodynamic stability derivatives from 
flight data. Several other methods are being used in system-identification problems; 
examples are given in references 3 and 4. 


Newton-Raphson Method 

The theoretical development of the Newton-Raphson method is based on a Taylor- 
series expansion of several variables. For example, the vector R(x) represents a set 
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of residue equations of the form 


R^(x) — Rj ^x^, X 2? • • • » x m 'j ~ ® 

R 2^ = R 2( X 1’ x 2’ ‘ * ' ’ x m) = 0 

r ( 14 ) 


R (x) = R (x . , x 0 , 
m ' mi 1 2 


m 


= 0 


J 


Expanding R(x) in a Taylor series about an approximate solution vector x n obtained 
in the nth iteration gives 



SR(x) 


3x J 


( x n+l - x n) + ■ • ' 


(15) 


If x „ - x is sufficiently small, the higher order terms of equation (15) can be 
n+l n 

neglected. Thus 



R(x) 


Jn 




(16) 


In the final iteration, the residue vector R(x) must become zero for the vector x to 
be a solution; equation (16) is then written as 



+ 



(17) 


8 



or as 




3R..(x) 

3R,(x) 

X 

3R,(x) 

X 




R x (x) 


X 




v _ 

Y 


3x^ 

3x 2 * * * 

3x 

m 


x l, n+1 

A 4 

l,n 

r 2 (x) 


3R 2 (x) 

3R 2 (x) 

BR 2 (x) 



x o ^ 


3 Xj 

3x 2 

3x 

m 


x 2, n+1 

• 

+ 

• 

• 

• 


• 


R (x) 
nr 1 


3R (x) 
m ' 

3R m (x) 

3R (x) 
m ' 


• 

Y 

X 

m, n 

n 


3x 2 

3 x 

m 


in, n+1 

- 

U 


n 

- 

- 


in matrix form. 

Denoting the Jacobian matrix by J 


3Rj(x) 

3Rj(x) 

BRjW 

3Xj~ 

3x 2 • • • 

3x 

m 

3R 2 (x) 

3R 2 (x) 

3R 2 (x) 

3Xj 

ax 2 

3x 

m 

3R (x) 
m x ' 

3R m® 

3R (x) 
m 

3Xj 

3x 2 • ’ ' 

3x 

m _ 


(18) 


(19) 


equation (18) is written in vector form as 


n+1 



( 20 ) 
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Equation (20) is the recurrence relation for the Newton-Raphson algorithm, which con- 
verges rapidly when the initial guess vector x is in the convergence space. Problems 
arise when irregularities occur in the individual equations. The Newton-Raphson meth- 
od has, for many problems, the properties of quadratic and monotone convergence 
(ref. 5). 

In an application of the Newton-Raphson algorithm, the nonlinear algebraic equa- 
tions (eq. (11)) are written as 


0 . 




R 1 (x) = y Q cos Xt 1 +-^~ sin Xt 1 - y x = 0 


— ■> o v 

R 2 (x) = y Q cos Xt 2 + -j- sin Xt 2 - y 2 = 0 / 


- y 0 
R 3 (x) = y Q cos Xt 3 +-^- sin Xt g - y 3 = 0 


( 21 ) 


The Jacobian matrix is 


[cos Xt|J [^ sin Xt[j j^y Q ti + -|^sin Xtj + — ^ cos Xt^ 
[cos Xt 2 ] [A sin Xt 2 J [- (y 0 t 2 + ~)sin Xt 


^ 0*2 . 

2 + x cos Xt 2 


[cos Xt 3 ] sin Xt 3 ] [-(y 0 t 3 +^) 


y 0*3 

Jsin Xtg + ^ cos X.t 3 


= J 


( 22 ) 


Using the recurrence equation (eq. (20)), x and R(x) are defined by 



y o 


I 

PS 

1 

X = 

y o 

R(x) = 

r 2 ( x ) 


X 


_R 3 (x)_ 


( 23 ) 
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Two initial guess vectors for x are selected to show the convergent properties of the 
problem. t Values of y^, y^, and X are presented in table I for each iteration of the 

initial estimate vectors of (0, 1, 1) and (0. 5, 0. 5, 0. 5). The following are boundary con- 


ditions to 

be satisfied: 


1 . 

y x = 1.00 

at 

t x = 0. 0 

2. 

y 2 = i - 36 

at 

t 2“ 0 * 5 

3. 

y 3 = 1.38 

at 

tg^.O 


TABLE I. - NEWTON- RAPHSON ITERATIVE SOLUTION TO EQUATIONS (11) 


State 

Initial 

Iteration 

variable 

guess 

1 

2 

3 

4 

5 

y o 

0. 0000 

1.0000 

-- 

-- 


-- 

y o 

1.0000 

1. 0000 

— 

-- 


-- 

X 

1.0000 

1.0000 

-- 

-- 

-- 

-- 

y o 

0. 5000 

0.98409 

0. 741428 

1.001728 

0.9999 

1.0000 

y o 

. 5000 

1.0000 

1. 00000 

1. 00000 

1. 0000 

1. 0000 

X 

. 5000 

1. 0000 

.92949 

1. 00861 

1. 0000 

1. 0000 


Newton- Raphson- Kantorovich Method 

The Newton- Raphson- Kantorovich method is an extension of the Newton-Raphson 
method to "function space" (ref. 6). The concept of function space, as explained by 
Lanczos (ref. 7), involves the replacement of a continuous function by a vector so that 
the vector describes the continuous function by a set of discrete points. 

The Newton- Raphson- Kantorovich method involves the solution of a set of linear 
differential equations with varying coefficients. In this set of equations, the solution of 
the linear differential equation converges, under appropriate conditions, to the solution 
of the nonlinear differential equations. Since the equations are linear and the principle 
of superposition applies, the boundary conditions can be satisfied at each iteration. 
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The development of the Newton- Raphson- Kantorovich algorithm starts with the 
consideration of the first-order vector differential equation 


Z = f(Z, t) 


( 24 ) 


where Z is composed of P dependent variables, and t is the independent time vari- 
able. Expanding equation (24) in function space gives 



n 


Neglecting the higher order terms and rewriting gives 



Equation (26) is now written in the form 


J n+1 


= A Z , 
n n+1 


+ B 


n 


where 


A 


n 




(26) 


(27) 


(28) 


The vectors Z and Z are linear with respect to the n + 1 solution. The A and 
B n terms in equation (27) are known functions that are calculated from the known nth 
solution or from the previous iteration. In the n+1 solution, the functions A r and 
B reflect the nonlinearity of the original differential equations. By using the recur- 
rence relation (eq. (27)), successive approximations are made to the nonlinear solu- 
tions until the desired accuracy is achieved. 
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Equation (27) is linear with varying coefficients and is easily solved through su- 
perposition. The particular solution is a solution of equation (27) with appropriate ini- 
tial conditions. The vector z|^(t) denotes this particular solution. The homogeneous 
solutions are governed by 

- v£i (t > lsisp < 29 > 


where P is the number of homogeneous solutions. The initial conditions for the homo- 
geneous solutions are selected so that the solutions will be linearly independent. The 
total solution of equation (27) is 


P 

5 n + l<‘> = * ZX 5 n*l (t) (3 ° 

i=l 

where the a. are specified to ensure satisfaction of the boundary conditions. The 
particular solution is Z^^(t), and the homogeneous solutions are the Z^j(t). 

The initial conditions for the particular and homogeneous solutions should show 
all the information known about the desired solution. Thus, the initial condition (t = 0) 
for the particular solution is 


z ( °) (t) = z (°\t) +'y%.z 

n+l v n w / j i n 
i=l 


(i) 


(t) 


(31) 


where the right-hand side of equation (31) is known from the previous iteration. The 
initial conditions for the homogeneous solutions are designated to be approximately the 
same, except that arbitrary perturbations of elements of these vectors are made to en- 
sure that 


det(zW( t ) 


5 n*X< l > 


) 7^0 

/t=o 


(32) 


These strategies ensure that — 0 as convergence is approached and therefore that 
a straightforward indication of convergence is obtained. 
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The number of boundary conditions is denoted by L. The individual boundary 
conditions are denoted by S k (1 < k < L). These boundary conditions are placed on 

elements of Z ^(t) and are satisfied by 


S k = [ Z n°i(‘k)] + S“i[ Z n’l(‘k)] 1 S k S L 

1 = 1 


( 33 ) 


where t, is the time at which S. is measured, 
k k 

If L = P, the following matrix equation yields the a. upon inversion of the coef- 
ficient matrix. 


Ca - d 

where 

C « = [ Z nil( h) 

defines the elements of the coefficient matrix. The elements of the right-hand vector 
are defined by 

d i - s i - 

If L > P and all boundary conditions are to be satisfied in a least-squares sense, the 
following matrix equation 


( 36 ) 


( 34 ) 


( 35 ) 


Ea = e 


( 37 ) 
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is used, where 


T 

E = C C 
— 

e = C d 


(38) 


When the appropriate matrix-inversion technique is used, the a's are easily deter- 
mined from equation (37). 


Once the vector a is determined, the solution and initial conditions for the non- 
linear differential equations are updated for the next iteration. The updated solution 
vector is calculated by using equation (31). The updated initial conditions are also cal- 
culated by using equation (31) and letting t = 0. The process is repeated with the up- 
dated solution and initial conditions until the desired accuracy is achieved. 

In an application of the Newton- Raphson- Kantorovich algorithm, the second-order 
differential equation is written as 


x + £x = 0 


(39) 


where the observed responses or boundary conditions are 


1 . 

x = Sj 

at 

t = t l 

2. 

x = s 2 

at 

t = t 2 

3. 

x = S 3 

at 

t ~ t 3 


The initial conditions Xq, Xq, and £ are now determined in order to satisfy the speci- 
fied boundary conditions. Equation (39) is reduced to two first-order differential equa- 
tions. 


x 


= u 


(40) 


and 


u = -4x 


(41) 
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Since | is a constant 


k = 0 


(42) 


Three first-order differential equations are then formed. 



X 


X 

Z 1~ Z 2 

z = 

u 

Z = 

u 

z 2 = -z 1 z 3 


k 


£ 

Z 3 = 0 


(43) 


Equation (39) is expanded according to equation (27). 


X n+1 U n+1 


(44) 


U n+1 ^n X n ^n( X n+l ” X n) ~ X n(^n+1 ^n) ~^n X n+l ^n+l X n + X n^n 


and 



(46) 


(Although equation (39) is linear, the same procedure applies if equation (39) is nonlin- 
ear. The procedure is independent of whether or not a closed-form solution exists. ) 
Before particular and homogeneous solutions can be determined, an initial guess must 

be made for the vector [Z(t)] t _Q. 


[z(t )] t=0 = 


U f 


(47) 
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Ill 



Time, l, sec 


Figure 1. - The first iterative solutions to 
equations (27) and (29) for the initial con- 
ditions defined in equations (48) and (49). 


An initial guess must also be made for the 
solution of the governing differential equa- 
tion (eq. (39)). A first approximation can 
be made by using linear interpolation be- 
tween the specified boundary conditions, as 
shown in figure 1. The assumed solutions 
and initial conditions are needed to calculate 
the n subscripted terms in equation (27). 

By using the initial conditions that re- 
sult in linearly independent solutions, the 
homogeneous and nonhomogeneous equations 
are numerically integrated over the time in- 
terval of interest such as t^ to t^. 

The initial conditions for the nonho- 
mogeneous solution (eq. (27)) are 



x o 


"1. 00000' 

Zj (0) (0) = 

u o 

= 

1. 00000 


J . 


_1.00000_ 


(48) 


The initial conditions for the homogeneous 
solution (eq. (29)) are 



0. 10000“ 


“1. 00000“ 


"1. 00000“ 

Zj (1) (0) = 

1. 00000 

Zj (2) (0) = 

0. 10000 

Zj (3) (0) = 

1. 00000 


_1. 00000 


_1. 00000_ 


_0. 10000 


(49) 


From equation (32), it is obvious that these initial conditions result in linear independ- 
ent solutions. The diagonal elements of the matrix of initial conditions for the homoge- 
neous solutions are perturbed by a factor of 0. 10000 for each iteration to ensure linear 
independence of the solutions. The perturbation factor is a specified parameter that 
can vary in magnitude depending on the type of problems being solved. The solutions 
to equations (27) and (29) for the previously discussed initial conditions are shown in 
figure 1. To ensure fast convergence, the selection of the starting values for the ini- 
tial conditions should reflect the best available values of the initial conditions. By 
using equation (33), three linear independent equations are formed. 
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§ 1 ■= CVl) - + + “s^nilOl) 

®2 “ 2 ^l( t 2) > (50) 

®3 = + “l^nVlt's) + a 2^l{h) + a ^nkh) ^ 

If the boundary conditions are specified displacements, the vector equation (50) is re- 
duced to a scalar equation. The scalars are the displacement component of the vector 

Z. In this instance, the following specified boundary conditions are displacements: 

1. Sj = 1. 00 at t 1 = 0. 00 

2. S 2 = 1. 00 at t 2 = 1. 00 

3. S 3 = -1.00 at t 3 = 1.50 

When the displacement component of the vector Z is used as specified in figure 1, the 
three linear independent equations (eq. (50)) are reduced to 


1. 00 = 1. 00 + 0. 10a x + 1. 00a 2 + 1. 00a 3 


1. 00 = 1. 38 + 0. 89a + 0. 62a 2 + 1. 89a 3 


1. 00 = 1. 04 + 1. Ola x + 0. 17a 2 + 2. 14a g 


Once the homogeneous and nonhomogeneous solutions have been determined, a standard 
matrix-inversion technique is used to determine the unknown a’s. After the a's have 
been determined, the solution to the governing differential equation is updated by using 
equation (31). The initial conditions and coefficients are also updated by using equa- 
tion (31). The process is repeated with the updated solution and initial conditions until 
convergence has been achieved. Values for x Q , Xq, and £ are presented in table II 

for each iteration on the solution to illustrate the convergent characteristics of this al- 
gorithm. 
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TABLE H. - NEWTON-RAPHSON-KANTOROVICH ITERATIVE 


SOLUTION TO EQUATION (39) 


State 

Exact 

Initial 



Iteration 





variable 

solution 

guess 

1 

2 

3 

4 

5 

6 

X 

o 

1. 00000 

1. 00000 

1. 00000 

1. 00000 

1. 00000 

1. 00000 

1. 

00000 

1. 

00000 

X 

0 

1. 20924 

1. 00000 

.684269 

1. 03947 

1. 19223 

1.20907 

1. 

20924 

1. 

20924 

1 

4. 38681 

1. 00000 

2.46593 

3. 85735 

4.33703 

4. 38633 

4. 

38681 

4. 

38681 


APPLICATIONS 


The quasilinearization techniques for solving system- identification problems can 
be applied to many categories of dynamics and related fields. The examples discussed 
in this report are confined to problems involving the dynamics of several one-degree- 
of-freedom systems. The numerical examples are designed to relate the accuracy of 
the observations of the pertinent state variable to the accuracy to which the coefficients 
and initial conditions of the governing differential equation are calculated. Specific ap- 
plications of a quasilinearization technique to the system- identification problems of 
parachute dynamics and automobile coasting dynamics are discussed. 


Numerical Examples 

The following nonlinear ordinary second-order differential equation is used for 
one set of numerical examples to relate the effect of noisy boundary conditions on the 
determination of the coefficients and initial conditions. 


H + A § + By + C&) 3 + Dy 3 = 0 (52) 

at \ / 


The independent variable t is time, and the dependent variable y is displace- 
ment. The initial conditions on y and dy/dt and the coefficients A, B, C, and D 
are selected so that oscillatory motion results. 
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If the coefficients C and D are set equal to zero, equation (52) reduces to the 
linear ordinary second-order differential equation 

^-J + A^ + By = 0 (53) 

dt 


which will be used as the governing differential equation for the second set of numerical 
examples. 

Equation (52) is analogous to a spring-mass system with linear and nonlinear 
properties of spring and damping forces. To develop precision data for the displace- 
ment history for equations (52) and (53), values for the coefficients and initial condi- 
tions were assigned and numerically integrated (ref. 8). The displacement data 
developed by numerical integration have six-digit accuracy. Noise was superimposed 
on the integrated displacement data by rounding off the six-digit-accurate displacement 
values to four- and two-digit accuracy. By making the numerically integrated data 
analogous to experimentally measured data, the Newton-Raphson-Kantorovich method 
is applied to the problem of identifying the differential- equation coefficients and initial 
conditions which best fit the experimental data. The degrees of accuracy of the numer- 
ically integrated data are assumed to be analogous to experimental error. 

Comparisons are presented in table m for the initial conditions and coefficients 
computed by using quasilinearization for the nonlinear differential equation (eq. (52)). 


TABLE m. - TABULATED VALUES OF INITIAL CONDITIONS 
AND COEFFICIENTS FOR EQUATION (52) a 


State variable 
or coefficient 

Exact 

solution 

6-digit 

accuracy 

4- digit 
accuracy 

2 -digit 
accuracy 


0. 00000 

-0. 000001 

0. 000010 

-0. 001527 

*0 

1. 00000 

1. 00002 

.999846 

1. 02469 

A 

. 10000 

. 100027 

. 100171 

. 059754 

B 

3. 00000 

3.00015 

2.99986 

2.99733 

C 

. 20000 

. 199912 

. 199583 

. 303364 

D 

4. 00000 

4. 00063 

4. 00348 

4.03972 


a Fourteen specified displacements obtained every 0. 5 second from t = 0 to 
t = 6. 5. 
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In this computation, six-, four-, and two-digit accuracy was used on the specified dis- 
placement history. The specified displacements were obtained at intervals of 0. 5 sec- 
ond from t=0 to t=6.5 seconds. For example, if the boundary conditions in 
in table 331 are specified to four-digit accuracy, the coefficient A is calculated to be 
0. 100171, as compared to the exact value of 0. 100000 for A. Similar comparisons for 
the linear differential equation (eq. (53)) are presented in table IV. 


TABLE IV. - TABULATED VALUES OF INITIAL CONDITIONS 
AND COEFFICIENTS FOR EQUATION (53) a 


State variable 
or coefficient 

Exact 

solution 

6 -digit 
accuracy 

4 -digit 
accuracy 

2-digit 

accuracy 

y 0 

0. 0000 

0. 000000 

0. 000000 

-0. 001870 

y 0 

1. 00000 

1.00002 

1. 00001 

1. 00012 

A 

. 100000 

. 099999 

. 100030 

. 102050 

B 

3. 00000 

3. 00015 

3. 00015 

3. 00104 


3. 

Fourteen specified displacements obtained every 0. 5 second from t = 0 to 


The data in tables III and IV indicate that the initial conditions are calculated to 
the same accuracy as the specified data for the examples studied. The coefficients, 
however, have the same accuracy as the specified data for the linear differential equa- 
tion (eq. (53)) and one less digit of accuracy for the nonlinear differential equation 
(eq. (52)). This difference in accuracy is probably attributed to the higher ratio of the 
specified displacement data to the differential- equation parameters (initial conditions 
and coefficients). The time span for which the displacement data are specified also in- 
fluences computed parameter accuracy. 


Parachute Dynamics 

The Newton- Raphson- Kantorovich method of solving system- identification prob- 
lems is also applied to parachute dynamics. Parachutes are used primarily as aero- 
dynamic deceleration devices. The normal sequence of parachute deployment is 
ejection from the payload, inflation or opening (which may be controlled or uncon- 
trolled), and subsequent descent. 

The parachute dynamics during inflation are difficult to model mathematically 
because of changing parachute geometry, virtual mass effects, et cetera. The descent 
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dynamics after the parachute has stabilized from the inflation process are more 
straightforward. The dynamics of the parachute trajectory are approximated with good 
accuracy by the second-order nonlinear differential equation 



+ C D S 



= g 


(54) 


where x is displacement (height) in feet, t is time in seconds, M is payload mass 
in slugs, C D S is a drag-area parameter in square feet, p is air density in slugs per 

2 2 

cubic feet, g is acceleration of gravity in feet per second per second, d x/dt is 
vehicle acceleration in feet per second per second, and dx/dt is vehicle velocity in 
feet per second. 

Experimentally measured values of position versus time for a free-falling, fully 
inflated parachute are presented in table V (ref. 9). These data are for the time span 
following the parachute inflation process. During this period, the parachute payload 
system decelerates from the opening velocity to the steady-state velocity, which occurs 
when aerodynamic drag is equal to parachute payload weight. The fall distance for 
which the parachute data are taken is approximately 1000 feet. Because of the small 
change in atmospheric density during this period, the density is assumed constant. 

With equation (54) as the governing differential equation, the coefficient and ini- 
tial conditions are determined by using quasilinearization so that the solution to the 
differential equation satisfies the experimental displacement data in a least-squares 
sense (minimizes the square of the deviation). The predicted values of parachute dis- 
placement, velocity, acceleration, and drag-area parameter CpS are presented in 

table V. The values selected for the drag-area parameter and initial conditions corre- 
late generally within 1 foot of the measured trajectory for a free-fall distance of ap- 
proximately 1000 feet. 

Several numerical experiments are presented that relate the accuracy of the cal- 
culated coefficient and the initial conditions of the differential equation to the accuracy 
of the specified boundary conditions (experimental data). The following differential 
equations and initial conditions are assumed and numerically integrated to develop pre- 
cise data for the boundary conditions. 


2 . 2 
iL* + 0. 24708 x l0- 3 (f) = 32. 17 


(55) 



443. 2 ft/sec 


(56) 
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= 0.0 


(57) 


x 


t=0 


Values of displacement and velocity are presented in table VI for the solution of equa- 
tion (55). The tabulated values have six-digit accuracy. 

TABLE V. - PARACHUTE DATA 3, 

C D S = 147. 78 ft 2 , p = 0. 1354 X 10 -2 slug/ft 3 , 

M = 405 slugs 


Time, 

sec 

Experimental 

displacement, 

ft 

Predicted 

displacement, 

ft 

Predicted 

velocity, 

ft/sec 

Predicted 

acceleration, 

ft/sec 2 

0 . 0 

0 . 0 

0 . 0 

443. 20 

-16. 35 

. 2 

87. 5 

88.3 

440. 00 

-15. 66 

.4 

176.0 

176.0 

436.94 

-14.99 

.6 

263.0 

263. 1 

434.00 

-14.36 

. 8 

349. 2 

349. 6 

431. 19 

-13. 76 

1.0 

434. 7 

435.6 

428. 50 

-13. 19 

1. 2 

520. 8 

521. 0 

425.91 

-12.64 

1.4 

606. 7 

606. 0 

423.43 

-12. 12 

1. 6 

691. 7 

690. 4 

421. 06 

-11. 62 

1.8 

775. 0 

774.4 

418. 78 

-11. 15 

2. 0 

857.6 

857.9 

416. 60 

-10. 70 

2. 2 

940. 2 

941. 0 

414. 50 

-10. 27 


a These data are a portion of an overall parachute payload trajectory in which 
several parachutes were deployed at different time sequences. These data were se- 
lected because the parachute payload dynamics can be approximated by equation (54). 
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TABLE VI. - SOLUTION TO EQUATION (55) 


Time, 

sec 

Displacement, 

ft 

Velocity, 

ft/sec 

0.0 

0. 000 

443. 200 

.2 

88. 317 

439.998 

.4 

176.008 

436.932 

.6 

263. 099 

433.996 

.8 

349. 615 

431. 183 

1.0 

435.580 

428.487 

1.2 

521. 017 

425.903 

1.4 

605.948 

423.426 

1.6 

690. 394 

421. 050 

1.8 

774.375 

418. 771 

2. 0 

857.909 

416. 584 

2.2 

941. 014 

414.486 


The numerically developed data in table VI were rounded off to simulate the ex- 
perimental error. The results of using data of ±0. 001-, ±0. 1-, and ±5-foot displace- 
ment or position accuracy to determine the accuracy of coefficients and initial 
conditions are presented in table VII. 

The data indicate, for the velocity and time span studied, that the coefficient of 
the velocity- squared expression and initial velocity can be calculated to within 1 per- 
cent for 5-foot position errors. By extrapolating this value to the parachute experi- 
mental data, the drag-area parameter C D S can be calculated to within 1 percent if the 

experimental position error is within 5 feet. 
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TABLE Vn. - EFFECT OF NOISY BOUNDARY CONDITIONS ON VALUES 
OF DIFFERENTIAL- EQUATION PARAMETERS 



Parameter 

Exact 

Specified boundary-condition accuracy, ft 

value 

±0. 001 

±0. 1 

±5.0 

o 

II 

443. 200 

443.200 

442.830 

443. 175 

c i 

0. 2470800 x 10" 3 

0. 2470818 x 10‘ 3 

0. 2453796 x 10‘ 3 

0. 2466957 x 10~ 3 


Coasting-Automobile Dynamics 

The dynamics of a coasting automobile are primarily influenced by two externally 
applied forces: aerodynamic drag and rolling friction. The aerodynamic drag is con- 
sidered to be a function of the square of the relative velocity of the automobile (relative 
to the wind). Rolling friction forces are assumed to be a function of the automobile 
weight. 

The governing differential equation for a coasting automobile for a no-wind con- 
dition is 


d 2 x 

dt 2 


pc 7 ,s 


^D°/dxV 
2M \dt) 


+ pg = 0 


(58) 


where p is the rolling friction coefficient. 

The differential equation for the coasting dynamics of an automobile is similar to 
the differential equation for a free-falling parachute. If p = -1, equation (58) is identi- 
cal to equation (54). The system- identification problem of the automobile includes the 
determination of rolling friction coefficient p. 

The experimental data for the coasting automobile are in the form of a velocity- 
time history (ref. 10). Table VIII presents experimentally measured values of velocity 
versus time for a coasting automobile. By using equation (58) as the governing differ- 
ential equation, the coefficients and initial conditions are determined such that a 
’’best" fit on the experimental data was obtained. Also presented in table VIII are the 
values of velocity, displacement, drag-area parameter C^S, and rolling friction coef- 
ficient p predicted by 'quasilinearization. To relate the accuracy of the predicted 
coefficients and initial conditions to the accuracy of the specified boundary conditions, 
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the differential equation and initial conditions 



1. 303065 x 10 



+ 0.6671862 = 0 


(59) 


x 


t=0 


0.0 ft 


(60) 



116.90158 ft/sec 


(61) 


are numerically integrated. The integrated values of displacement and velocity are 
considered to be analogous with experimental data of six-digit accuracy. 

TABLE VIE. - COASTING-AUTOMOBILE DATA 

x t=Q = 0, (dx/dt) t=Q = 116. 90158 ft/sec, C D S = 8. 80719 ft 2 , 
p = 0. 00238 slug/ft 3 , p = 0. 0207201, M = 80. 43 slugs . 


Time, 

sec 

Experimental measured 
velocity, ft/sec 

Predicted 
velocity, ft/sec 
(a) 

Predicted 
displacement, ft 
(a) 

0 

117.3 

116.902 

0. 000 

5 

104.9 

105. 520 

555. 369 

10 

95. 3 

95. 604 

1057.641 

15 

87.3 

86.852 

1513. 350 

20 

79.2 

79.038 

1927. 725 

25 

71.9 

71.990 

2305.007 

30 

65. 7 

65.574 

2648.678 

35 

59.8 

59.684 

2961.622 

40 

54.0 

54. 235 

3246. 250 


a These numerically integrated values have six-digit accuracy. 
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Numerical experiments were conducted to determine the influence of experimental 
error on the predicted differential- equation coefficient and initial conditions. Two types 
of numerical experiments were conducted: in one type of experiment, the velocity data 
were used as the specified boundary condition, and in the other type of experiment, the 
displacement data were used as the specified boundary conditions. These experiments 
were designed to compare the accuracy in computing differential-equation coefficients 
and initial conditions for specified boundary conditions of different orders. The exper- 
imental error was developed by rounding off the integrated data to accuracies of 1 and 
5 ft/sec for velocity and to accuracies of 1, 5, and 10 feet for displacement. 

Values of initial conditions and coefficients are presented in table IX for each set 
of data at a given experimental error. For a ±5 -foot displacement error, the co eff i- 

_4 

cient C is computed to be 1. 30469 x 10 , which is within 0. 12 percent of the exact 

value of CL. 

An error of ±10 feet in experimental displacement data results in approximately 
the same accuracy in computation of the differential- equation parameters as a ±5-ft/sec 
error in the experimental velocity data. An error of ±5 feet in the displacement data 
gives better accuracy in the computed differential- equation parameters than a ±l-ft/sec 
error in the velocity data. 
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TABLE IX. - EFFECTS OF NOISY BOUNDARY CONDITIONS ON CALCULATED PARAMETERS OF DIFFERENTIAL EQUATION (58) 


^T +C l(ar) + C 2 = °’ P = 0.00238 slug/it 3 . CpS = 2CjM/p. 
M = 80.43 slugs. H - C^g 


Parameter 

Exact value 

Specified boundary-condition accuracy 

f Displacement ; 

Velocity 

±1 ft 

Percent 

error 

(a) 

±5 ft 

Percent 

error 

(a) 

±10 ft 

Percent 

error 

(a) 

±1 ft sec 

! 

Percent 

error 

(a) 

±5 ft sec 


Percent 

error 

(a) 

(dx/dt)^. 

116.902 

116.840 

0.053 

117. 147 

0.21 

117.402 

0.42 

117. 124 


0. 19 

115.534 


1. 18 

ft sec 














c i 

1.30307 x 10' 4 

1.28490 x io' 4 

1.41 

1.30469 x 10~ 4 

0. 12 

1. 3941 1 x 10 -4 

6.53 

1. 27184 v 

io - 4 

2.45 

1.34169 ' 

io - 4 

5.69 

C 2 

0.667186 

0.679847 

1. 86 

0. 678837 

1. 72 

0.608339 

9.67 

0. 690197 


3. 33 

0. 586287 


13.80 

c D s, ft 2 

0.509086 

0.50199 

1. 41 

0.509723 

0. 12 

0. 544657 

6.53 

0. 496900 


2. 45 

0.539806 


5.69 

M 

i 

2.07201 x 10" 2 

2. 11133 x 10" 2 
.1 

1. 86 

2. 10819 x 10‘ 2 

1. 72 

1.88925 ■ 10' 2 

9.67 

1 

2. 14347 ■ 

,i 

io -2 

3. 33 

1. 82077 

io - 2 

13.80 


a 


Percent error = 


exact value 
computed value 


x 10‘ 2 


CONCLUDING REMARKS 


The quasilinearization techniques provide a systematic iterative approach to fit- 
ting a governing differential equation (or equations) to a set of measured observations. 
The advantage of this technique in system identification is that it permits the selection 
of the state variables that describe the measured observation whether the observation 
is displacement, velocity, or any other variable. The measured observations can be 
made from the standpoint of economy, accuracy, or convenience. In terms of economy, 
this technique may reduce instrumentation requirements. 

The accuracy of the predicted differential- equation coefficients and initial condi- 
tions is a function of experimental error -in the measured observations. Because of the 
nature of measuring devices, the higher the order of the state variable (such as velocity 
and acceleration) being measured, the noisier the experimental data. The quasi- 
linearization technique allows the selection of the state variable with the least noise or 
the lowest order, which is displacement in dynamics problems. 

The quasilinearization techniques are applied to the identification of several one- 
degree-of -freedom systems. Linear and nonlinear systems are used to relate the ac- 
curacy of the differential- equation coefficients and initial conditions for a given 
experimental error in boundary conditions. The data indicate that the initial conditions 
can be calculated to the same accuracy as the specified boundary conditions (e. g. , two- 
digit-accurate boundary conditions for two-digit-accurate initial conditions), whereas 
the computed coefficients have less accuracy than the specified boundary conditions. 

Specific system-identification problems of the free-falling parachute are dis- 
cussed. Results indicate, for the conditions studied, that the parachute aerodynamic 
drag-area parameter can be calculated with good accuracy during a transient, free- 
falling condition. The numerical examples indicate that the aerodynamic drag-area 
parameter can be calculated to ±1 percent when errors in parachute position data are 
as large as 5 feet for the parachute trajectory studied. 

Experimental data on the coasting automobile are presented in the form of a ve- 
locity history. The system-identification problem for a coasting automobile includes 
determination of the initial conditions, aerodynamic drag-area parameter, and rolling 
friction coefficient. Numerical experiments were conducted to relate errors in exper- 
imental velocity histories and displacement histories to the predicted accuracy of the 
differential- equation parameters. The results indicate that determination of the 
differential- equation parameters is almost twice as accurate when displacement data 
rather than velocity data with the same order of accuracy are used. 


National Aeronautics and Space Administration 
Manned Spacecraft Center 

Houston, Texas, April 25, 1969 
914-13-20-17-72 
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